####################################################################
# AGER FIRE GENERATOR
####################################################################
# Created: 11.1.2011
# Latest Updates: 4.2.2013 (Update: replaced wind direction code)
####################################################################
    ##METHODS:
          # Read MC1 monthly ERC values
          # Generate daily ERC values from monthly values using LOG std dev using B. Johnson methods
          # Determine the probability of a fire from Preisler smooth surface fit and logit-line
          # Determine the size of the fire from Greenfell splines and exp distribution
          # Determine wind speed and wind direction from Wind Rose Gust data - corrected!
          # Convert the fire size to burn period based Ager calibration data
          # Determine fuel moisture values from regression, and generate fuel moisture files
         ## Implement agricultural/grassland modifications to fuel moisture files ##
          # Create full fire list 
          # Create fire list file with fires > 1 ha
          # Create fire list file with fires > 1 ha AND implement ignition probability switch

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#  INFORMATION YOU NEED TO RUN PROGRAM:
#   1) Set filename path (working directory); place all files in this location (line 68)
#   2) Install necessary packages (lines 63-65) and load libraries (lines 71-73)
#   3) Files: a) monthly MC1 input data file (vcfutercgpred.csv)
#             b) historical fire events file (IgnitionProbData.csv)
#             c) historical fire ignitions file (erc-avgfiresize.csv)
#             d) hitorical wind gust file (BrushCreek_WindRose.csv)
#             e) hitorical wind direction file (BrushCreek-WindRose-Summary.csv)
#             f) calibration burn period and windspeed data (new_calibration_firesize_BP_july_2012.csv)
#             g) village creek historical fuel moisture file (villagecrk_fuelmoisture.csv)
#             h) oak knoll historical fuel moisture file (oakknoll_fuelmoisture.csv)
#   4) Start and end dates for future MC1 monthly data
#   5) Size of the study area, line 219, Eugene = 821 sq km
#   6) Set location to output fuel moisture files, see line 647
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#   CHOOSE YOUR CLIMATE SCENARIO (**search on keyword 'climate'**)
#---------------------------------------------------------------------
# You must specify the climate scenario in the following locations:
#   1) Line 86: select model by specifying column number
#   2) Line 127: select model by specifying column number
#   3) Line 647: specify path for fuel moisture files for Envision
#   4) Line 1604: change the firelist filename based on model for full list
#   5) Line 1622: change the firelist filename based on model with small fires removed
#   6) Line 1645: change the firelist ignition switch filename based on model
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#############################################
#  CHOOSE CRAN MIRROR
#--------------------------------------------
# Either select the CRAN mirror manually via your GUI interface or,
# run code below!

chooseCRANmirror(graphics=getOption("menu.graphics"))   ## RUN ONLY ONCE ##

######################################
#  PACKAGES NEEDED:
######################################
#   1. DATA.TABLE
######################################

install.packages("data.table")   ## RUN ONLY ONCE ##
install.packages("lattice")      ## RUN ONLY ONCE ##
install.package("mgcv")          ## RUN ONLY ONCE ##

#Set working directory
setwd("J:\\runfiregen\\")

#Load libraries
library (data.table)
library (lattice)
library (mgcv)

##################################################
# GENERATE DAILY MC1 ERC VALUES FROM MONTHLY MC1 ##############################
##################################################
  #To generate daily values within R....
  #Generate list of all years, months and days for our data record
  #For each day calculate the variance in ERC

mc1.month.mat <- read.csv("vcfutercgpred.csv") #Read monthly ERC values

## ---> Specify climate model!

erc.pred.mc1 <- mc1.month.mat[,4] #Define ERC variable where col4 = HADLEY, col5 = CSIRO, col6 = MIROC
yr.mc1 <- mc1.month.mat[,3] #Define year
mo.mc1 <- mc1.month.mat[,2] #Define month

#------------------------------------------
#Create data set with all years, months and days......

#Enter start and end dates for sequence
dateseq <- seq(as.Date("2007/1/1"), as.Date("2099/12/31"), by="day")

jday <- strftime(dateseq, format="%j") #Create julian day-of-year variable
month <- strftime(dateseq, format="%m") #Define month variable
month <- as.integer(month)
day <- strftime(dateseq, format="%d") #Define day variable
day <- as.integer(day)
year <- strftime(dateseq, format="%Y") #Define year variable (with century)
year <- as.integer(year)
fireday <- strftime(dateseq, format="%B-%d") #Define fireday variable

final.dates <- cbind(dateseq, year, month, day, jday)
final.dates <- as.data.frame(final.dates)
#---------------------------------------------

#Merge final dates file with original monthly ERC file
mc1.month.mat2 <- mc1.month.mat[,-1] #First delete id variable
mc1.daily.mat <- merge(final.dates, mc1.month.mat2, by=c("year", "month"), all=TRUE)
mc1.obs <- nrow(mc1.daily.mat) #Get number of observations
#There are 33,968 observations in the dataset

#R converts year and month to factor during the merge step- these must be converted back to maintain the correct
#vector order
mc1.daily.mat$year.int <- as.numeric(levels(mc1.daily.mat$year)[mc1.daily.mat$year]) #Convert factor to integer!
mc1.daily.mat$mon.int <- as.numeric(levels(mc1.daily.mat$month)[mc1.daily.mat$month]) #Convert factor to integer!
mc1.daily.mat.final <- mc1.daily.mat[,-1] #Remove old year factor variable
mc1.daily.mat.final <- mc1.daily.mat.final[,-1] #Remove old month factor variable

#Resort matrix to original order by year, month and day
mc1.daily.sort <- mc1.daily.mat.final[order(mc1.daily.mat.final$year.int, mc1.daily.mat.final$mon.int, mc1.daily.mat.final$day) , ]

## ---> Specify climate model!

erc.pred.mc1 <- mc1.daily.sort[,4] #Define ERC variable again where col4 = HADLEY, col5 = CSIRO, col6 = MIROC

#Calculate standard deviation of ERC (Equation from B. Johnson)
erc.sd <- -1.80293+3.22763*log(erc.pred.mc1+0.1) #OLD METHOD 
#erc.sd <- (-0.0099*(erc.pred.mc1*erc.pred.mc1)) + (0.6446*erc.pred.mc1) - 0.1824 ## NEW ##
#erc.sd <- ifelse(erc.sd <8, 8, erc.sd) ## NEW ##
  
#Calculate daily ERC (Method from B. Johnson)
#set.seed(7) #Inactivated to develop variability and multiple list files
erc.daily <- erc.pred.mc1 + (erc.sd*rnorm(mc1.obs)) #Specify # of observations

#Replace negative values with zero and round
erc.mc1 <- round((ifelse (erc.daily<0, 0, erc.daily)), digits=4)

#Create month vector for Wind Rose
month.sim <- mc1.daily.sort[,8]

##NOTE: Final ERC vector = erc.mc1

############################
# DETERIMINE FIRE IGNITION ###################################
############################
  #Based on ERC for 4 RAWS stations and FPA Ignition data for ignition study area, plus
    #Random sample of grid points for Ignition Training Area
  #RAWS Stations (North to South):
        #Village Creek
        #Mt. Yoncalla
        #Silver Butte
        #Provolt Seed Orchard
  #There are 2,184 missing dates between 8/5/1988 and 12/31/2008, thus only 5,270 days in record
  #No fires occured on missing RAWS dates
  #Data set contains ERC values for a given RAWS Station Zone matched spatially to ignitions
  #Thus dates may repeat for days with fires in multiple stations
  #For NON Ignitions three sets of ERC values were randomly selected for each day from grid
  #Data set contains date, RAWS Station (or region), historical daily ERC value, regional mean ERC,
      #fire event, year, month and julian day of year
#------------------------------------------------------------------------
#   STUDY AREA INFORMATION
#------------------------------------------------------------------------
#   Ignition Area: 
#           3,879,338 hectares 
#           4,222 ignitions, 201 ignitions/year
#           Average # of days per/year with >=1 ignition = 66 (35-160)
#   Lebanon Area: 
#           101,887 hectares 
#           2.6% of ignition area
#           45 ignitions, 2.1 ignitions/year
#           Fraction of observed fires in entire ignitions area = 0.012
#           Mean historical ERC = 34 (for ignitions)
#   Eugene Area: 
#           82,100 hectares 
#           2.1% of ignition area 
#           56 ignitions, 2.6 ignitions/year
#           Fraction of observed fires in entire ignition area = 0.013
#           Mean historical ERC = 32 (for ignitions)
#   Village Creek Historical Mean ERC = 27 (excluding zeros)
#   Mt. Yoncalla Historical Mean ERC = 24 (excluding zeros)
#   Silver Butte Historical Mean ERC = 27 (excluding zeros)
#   Provolt Seed Orchard Historical Mean ERC = 29 (excluding zeros)
#------------------------------------------------------------------------

erc.mat <- read.csv("IgnitionProbData.csv") #Historical ignition data and random sample of non-ignitions

erc.hist <- erc.mat[,2] #Define historical ERC variable
site <- factor(erc.mat[,3]) #Define RAWS station zones
resp <- erc.mat[,4] #Define response, event 0 = no fire, 1 = fire
ig.year <- erc.mat[,5] #Define ignition year
ig.month <- erc.mat[,6] #Define ignition month
ig.jday <- erc.mat[,7] #Define ignition day of year

fit1 <- bam(resp~te(erc.hist,ig.jday, k=c(3,5)), binomial) #Fit a smooth surface 

pi <- 5270*3/(46701*5270) #Calculate proportion sampled; 46,701 potential grid points in the ignition training area

#.......................................................................................................
#NOTE: Grid points from the Ignition training area were randomly selected for each day in the historical
#date record (5,270 days), with replacement, 3 times. These points were then matched to the RAWS zone
#in which they fell to determine the ERC value on a given day.  These represent NON-Ignitions
#.......................................................................................................

jday <- as.integer(jday) #Convert day of year from string to integer

nerc <- data.frame(erc.hist=erc.mc1, ig.jday=jday) #Where erc.mc1 = predicted future ERC
                                                   #Where jday = day of year for future dates

lp <- predict(fit1, nerc)-log(1/pi) #linear predictor (logit-line) for new erc values   
                                    #adjusted for sampling proportion

ep <- exp(lp)/(1+exp(lp))   #predicted probability of an ignition for the given erc values
                            #these probabilities are per day per grid (km^2)

EP <- 821*ep #predicted number of ignitions for the eugene study area, km^2

plot (erc.mc1, EP, main = "Predicted Ignition Probability - All Data",
      ylab = "Predicted number of ignitions")

ig.dat <- as.data.frame(cbind(jday, erc.mc1, EP)) #Create dataframe for graphing
jday200 <- as.data.frame(ig.dat[ig.dat$jday == 200, ]) #Create dataset for jday 200
jday300 <- as.data.frame(ig.dat[ig.dat$jday == 300, ]) #Create dataset for jday 300

#Plot data to compare how day of year changes ignition probability
plot (jday200$erc.mc1, jday200$EP, main = "Predicted Ignition Probability - JDay=200",
      ylab = "Predicted number of ignitions")
plot (jday300$erc.mc1, jday300$EP, main = "Predicted Ignition Probability - JDay=300",
      ylab = "Predicted number of ignitions")

# Equations created by Haiganoush Preisler 8.22.12
#-------------------------------------------------------------------

pred.fire.mc1 <- as.numeric(round(EP, digits=9)) #Redefine probability variable

#fire.prob <- cbind(erc.mc1, pred.fire.mc1) #Create matrix of mc1 data and fire probs

  ##NOTE: Final fire probability vector = pred.fire.mc1
  ##NOTE: Final fire probability matrix = fire.prob

########################
# DETERIMINE FIRE SIZE ############################
########################
#Code written by Isaac Greenfell
#Establish relationship between ERC and Fire Size from Historical Data

hist.mat <- read.csv("erc-avgfiresize.csv") #Historical FPA ignition data

ercseq <- hist.mat[,3] #Define ERC variable
fsizeseq <- hist.mat[,7] #Define Fire Size variable

#Convert fire size in acres to hectares
fsizeseq <- fsizeseq * 0.404685642

par(mfrow=c(1,1))

## SMOOTH SPLINE - NonParametric Regression Model
#Changing df will adjust how wiggly the line is
#Relationship between Fire Size and ERC
sm1 <- smooth.spline(fsizeseq~ercseq, df = 5)

predseq <- 1:120 #Domain over which fit is to be made

pred.sizeseq <- predict(sm1, predseq) #Calculate fitted values over range of ERC

plot(ercseq, fsizeseq)
lines(pred.sizeseq$x, pred.sizeseq$y, col="green")

#Assuming there is a fire using MC1 derived daily ERCs... 
mean.simseq <- predict(sm1, erc.mc1) #Calculate fitted values over range of MC1 ERC values
plot(mean.simseq$x, mean.simseq$y, 
     xlab="ERC", ylab="fire size (ha)") #simseq$x = ERC, simseq$y = firesize

#Generate random number from exponential distribution with rate
#This is conditional on day
out.sim <- rexp(erc.mc1, 1/mean.simseq$y)

#Plot new simulation data versus historical
par(mfrow = c(2,1))

max.ha <- 800

plot(mean.simseq$x, out.sim, ylim = c(0, max.ha))
plot(ercseq, fsizeseq, ylim = c(0, max.ha))

fire.size <- round(out.sim, digits=2) #Redefine variable

##NOTE: Final fire size vector (ha) = fire.size

##########################################################
# GENERATE WIND SPEED AND AZIMUTH FROM WIND ROSE GUST   #####################
##########################################################
#Based on Brush Creek RAWS Station historical data
#Years 2000 to 2010, units = mph
#Data are daily summary time series or wind rose summaries

wind.rose <- read.csv("BrushCreek_WindRose.csv") #Read wind data
rose.obs <- nrow(wind.rose)

w.date <- as.Date((wind.rose[,1]), "%m/%d/%Y") #Define date
w.month <- strftime(w.date, format="%m") #Define month variable
w.year <- as.integer(strftime(w.date, format="%Y")) #Define year variable (with century)
w.myr <- strftime(w.date, format="%m-%Y") #Create Month-by-Year variable

wind.rose$w.myr <- w.myr #Add month-by-year variable to dataframe

w.rose <- as.data.frame(wind.rose[, c(7,8)]) #Select needed variables

w.rose.mean <- aggregate(w.rose, list(as.factor(w.rose$w.myr)), mean, na.rm = T) #Calculate mean wind gust by month-year
w.rose.mean$w.mo <- as.integer(substr(w.rose.mean$Group.1, 1, 2)) #Create month variable

w.rose.sm <- as.data.frame(w.rose.mean[, c(2, 4)]) #Create dataframe with only needed variables; NULLs CAUSE WARNINGS
w.rose.sm[10,1] <- 24 #Replace missing values for one year with average of other years -January
w.rose.sm[21,1] <- 26.68 #Replace missing velus for one year with average of other years -February

#Create separate arrays for each month
w.jan <- subset(w.rose.sm, w.mo == 1)
w.jan <- as.vector(w.jan[,1])
w.feb <- subset(w.rose.sm, w.mo == 2)
w.feb <- as.vector(w.feb[,1])
w.mar <- subset(w.rose.sm, w.mo == 3)
w.mar <- as.vector(w.mar[,1])
w.apr <- subset(w.rose.sm, w.mo == 4)
w.apr <- as.vector(w.apr[,1])
w.may <- subset(w.rose.sm, w.mo == 5)
w.may <- as.vector(w.may[,1])
w.jun <- subset(w.rose.sm, w.mo == 6)
w.jun <- as.vector(w.jun[,1])
w.jul <- subset(w.rose.sm, w.mo == 7)
w.jul <- as.vector(w.jul[,1])
w.aug <- subset(w.rose.sm, w.mo == 8)
w.aug <- as.vector(w.aug[,1])
w.sep <- subset(w.rose.sm, w.mo == 9)
w.sep <- as.vector(w.sep[,1])
w.oct <- subset(w.rose.sm, w.mo == 10)
w.oct <- as.vector(w.oct[,1])
w.nov <- subset(w.rose.sm, w.mo == 11)
w.nov <- as.vector(w.nov[,1])
w.dec <- subset(w.rose.sm, w.mo == 12)
w.dec <- as.vector(w.dec[,1])

gust <- rep(1, mc1.obs) #Define variable

#Loop through each observation to randomly select max. wind gust by month
for (i in 1:mc1.obs){
  
  if (month[i] == 1) {
    gust[i] <- sample(w.jan, size=1);
  } else if (month[i] == 2) {
    gust[i] <- sample(w.feb, size=1);
  } else if (month[i] == 3) {
    gust[i] <- sample(w.mar, size=1);
  } else if (month[i] == 4) {
    gust[i] <- sample(w.apr, size=1);
  } else if (month[i] == 5) {
    gust[i] <- sample(w.may, size=1);
  } else if (month[i] == 6) {
    gust[i] <- sample(w.jun, size=1);
  } else if (month[i] == 7) {
    gust[i] <- sample(w.jul, size=1);
  } else if (month[i] == 8) {
    gust[i] <- sample(w.aug, size=1);
  } else if (month[i] == 9) {
    gust[i] <- sample(w.sep, size=1);
  } else if (month[i] == 10) {
    gust[i] <- sample(w.oct, size=1);
  } else if (month[i] == 11) {
    gust[i] <- sample(w.nov, size=1);
  } else if (month[i] == 12) {
    gust[i] <- sample(w.dec, size=1)
  }
} # ends loop

w.spd <- round(gust, digits=0) #Redefine wind gust as w.spd

#-------------------------------------------------------------------
#Wind Direction Data: Summary of dominant wind direction by month (percent)

wr.df <- read.csv("BrushCreek-WindRose-Summary.csv", row.name = 1) #Read wind rose summaries by month

wr.df <- wr.df/100 #Calculate probability from percent

##Note: That a calm threshold was set when windrose data were generated.
#       Calm threshold = 2.7 mph. Probabilities need to be corrected to sum to 1.

mo.prob <- apply(wr.df, 1, "sum") #Calculates sum of probabilities across months

wr.corr <- wr.df/mo.prob #Correct probabilities with sum by month

check <- apply (wr.corr, 1, "sum") #Check that sum of probabilities are now close to 1

sim.mo <- as.character(month.sim) #Convert month to character needed by for loop

dir.sample <- colnames(wr.df) #Grab direction row names

dir.sim <- rep(NA, mc1.obs ) #Create array of 1000 dummy samples for simulated direction

for(i in 1:mc1.obs)
{
  cur.mo <- sim.mo[i] #set cur.mo to vector of months for simulation period
  
  row.prob <- wr.corr[cur.mo,] #Set row probability to matching prob in the matrix
  
  row.sum <- sum(row.prob) #Calculate the sum of the probs of that row (should be 1)
  
  if(row.sum > 0)
  {  
    row.prob <- row.prob
    
    dir.sim[i] <- sample(dir.sample, 1, prob = row.prob)
    
  }
}   # ends loop


dir.substr <- substr(dir.sim, 2, nchar(dir.sim)) #Take direction value from string
azimuth <-as.integer(dir.substr) #Convert direction to integer

##NOTE: Final wind speed vector = w.spd
##NOTE: Final wind direction vector = azimuth

##############################
# RUN BURN PERIOD TRANSLATOR #####################
##############################
# Calibration data created by A. Ager, 7.25.12

#Data set with burn period and fire size (ha)

calibrate <- read.csv("new_calibration_firesize_BP_july_2012.csv") #Read calibration data

bp <- calibrate[,1] #Define burn period variable
Hectares <- calibrate[,2] #Define fire size (hectares) variable

par(mfrow = c(1,1)) #Set graph output parameters
plot(bp ~ Hectares) #Graph relationship

bp.reg <- lm((bp~Hectares+I(Hectares^2)-1)) #Polynomial linear regression model forcing intercept through 0,0
summary(bp.reg)

sizeha <- fire.size #Redefine fire.size data

preddat <- data.frame(Hectares = sizeha) #Create new dataframe with redefined fire size

pred.bp <- predict(bp.reg, newdata = preddat) #Predict new burn period with calibration data

plot(preddat[,1], pred.bp,
     xlab= "predicted fire size (ha)",
     ylab= "predicted burn period (min)") #Plot fire size by predicted burn period

burnp <- round(pred.bp, digits=0) #Redefine burnperiod variable and round

##NOTE: Final burn period vector = burnp

####################################
# GENERATE FUEL MOISTURE FILE LINK ##########################
####################################
#Village Creek ERC-fuel moisture data from Charles McHugh
    # villagecreek_frisk_file_May_Oct.txt
#Oak Knoll ERC-fuel moisture data from Charles McHugh
    # oak_knoll_frisk_file_may_oct.txt
#IMPORTANT NOTE: In order to capture high historical ERCs, Oak Knoll data will be used to predict
  # ERC Values above 60!!

fuel <- read.csv('villagecrk_fuelmoisture.csv')#Village Creek Fuel Moisture

erc.fuel <- fuel[,1] #Define ERC variable
fm1 <- fuel[,2] #Define 1 hr fuel moisture variable
fm10 <- fuel[,3] #Define 10 hr fuel moisture variable
fm100 <- fuel[,4] #Define 100 hr fuel moisture variable
fm1000 <- fuel[,5] #Define 1000 hr fuel moisture variable
fmherb <- fuel[,6] #Define herbaceous fuel moisture variable
fmwood <- fuel[,7] #Define woody fuel moisture variable

fuel.oak <- read.csv('oakknoll_fuelmoisture.csv')#Oak Knoll Fuel Moisture
fuel.high <- fuel.oak[fuel.oak$ERC>=60, ]#Grab only data where ERC >=60

erc.fuel.oak <- fuel.high[,1] #Define ERC variable -- Oak Knoll
fm1.oak <- fuel.high[,2] #Define 1 hr fuel moisture variable -- Oak Knoll
fm10.oak <- fuel.high[,3] #Define 10 hr fuel moisture variable -- Oak Knoll
fm100.oak <- fuel.high[,4] #Define 100 hr fuel moisture variable -- Oak Knoll
fm1000.oak <- fuel.high[,5] #Define 1000 hr fuel moisture variable -- Oak Knoll
fmherb.oak <- fuel.high[,6] #Define herbaceous fuel moisture variable -- Oak Knoll
fmwood.oak <- fuel.high[,7] #Define woody fuel moisture variable -- Oak Knoll

## FUEL MOISTURE - 1 HR FUELS------------------------
fm1.reg <- lm((fm1~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fm1.reg)

fm1.reg.oak <- lm((fm1.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model-- Oak Knoll
summary(fm1.reg.oak)

par(mfrow = c(4,1)) #Set graph output parameters
max.erc <- 120 #Set x-axis maximum value

plot(erc.fuel, fm1, xlim = c(0, max.erc), main = "Village Creek 1hr Fuel Moisture") #Plot Village Creek
lines(erc.fuel, predict(fm1.reg), col="blue") 

plot(erc.fuel.oak, fm1.oak, xlim = c(60, max.erc), main = "Oak Knoll Fuel Moisture- High ERC") #Plot Oak Knoll
lines(erc.fuel.oak, predict(fm1.reg.oak), col="blue")

tempx <- erc.mc1 #Redefine erc.mc1 data
hightempx <- subset(erc.mc1, erc.mc1 > 59.9999) #Redefine erc.mc1 data -- Oak Knoll

newdat <- data.frame(erc.fuel = tempx) #Create new dataframe with redefined mc1
newdat.high <- data.frame(erc.fuel.oak = hightempx) #Create new dataframe with redefined mc1 -- Oak Knoll

pred.fm1 <- predict(fm1.reg, newdata = newdat) #Predict new 1hr fuel moistures with erc.mc1 data
pred.fm1.oak <- predict(fm1.reg.oak, newdata=newdat.high) #Predict new 1hr fuel moistures with erc.mc1 data --  Oak Knoll

plot(newdat[,1], pred.fm1, main = "Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 1hr fuel moisture
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fm1.oak, main = "Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 1hr fuel moisture

## FUEL MOISTURE - 10 HR FUELS------------------------------------------------------------------------
fm10.reg <- lm((fm10~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fm10.reg)

fm10.reg.oak <- lm((fm10.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model
summary(fm10.reg.oak)

plot(erc.fuel, fm10, xlim = c(0, max.erc), main = "Village Creek 10hr Fuel Moisture")#Plot Village Creek
lines(erc.fuel, predict(fm10.reg), col="blue")

plot(erc.fuel.oak, fm10.oak, xlim = c(60, max.erc), main = "Oak Knoll Fuel Moisture- High ERC")#Plot Village Creek
lines(erc.fuel.oak, predict(fm10.reg.oak), col="blue")

pred.fm10 <- predict(fm10.reg, newdata = newdat) #Predict new 10hr fuel moistures with erc.mc1 data
pred.fm10.oak <- predict(fm10.reg.oak, newdata=newdat.high)#Predict new 10hr fuel moisture with erc.mc1 data -- Oak Knoll

plot(newdat[,1], pred.fm10, main = "Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 10hr fuel moisture
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fm10.oak, main = "Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 10hr fuel moisture - Oak Knoll

#Fuel moistures above 2.7 are incorrectly predicted!  These data will be replaced.
#Fuel moistures should reach asymptote and not continue to drop.
#Replace values where fuel moisture <2.7 with 2.69 based on historical Oak Knoll data

pred.fm10.oak2 <- ifelse(pred.fm10.oak <2.7, 2.69, pred.fm10.oak) #Correct asymptote

plot(newdat.high[,1], xlim = c(60, max.erc), pred.fm10.oak2, main = "Predicted Fuel Moisture for High ERC - corrected",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 10hr fuel moisture - Oak Knoll

## FUEL MOISTURE - 100 HR FUELS------------------------------------------------------------------------
fm100.reg <- lm((fm100~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fm100.reg)

fm100.reg.oak <- lm((fm100.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model-- Oak Knoll
summary(fm100.reg.oak)

plot(erc.fuel, fm100, xlim = c(0, max.erc), main="Village Creek 100hr Fuel Moisture")#Plot Vilage Creek
lines(erc.fuel, predict(fm100.reg), col="blue")

plot(erc.fuel.oak, fm100.oak, xlim = c(60, max.erc), main="Oak Knoll Fuel Moisture- High ERC")#Plot Vilage Creek
lines(erc.fuel.oak, predict(fm100.reg.oak), col="blue")

pred.fm100 <- predict(fm100.reg, newdata = newdat) #Predict new 100hr fuel moistures with erc.mc1 data
pred.fm100.oak <- predict(fm100.reg.oak, newdata = newdat.high) #Predict new 100hr fuel moistures with erc.mc1 data-- Oak Knoll

plot(newdat[,1], pred.fm100, main="Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 100hr fuel moisture
plot(newdat.high[,1], pred.fm100.oak, xlim = c(60, max.erc), main="Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted 100hr fuel moisture -- Oak Knoll

## FUEL MOISTURE - 1000 HR FUELS------------------------------------------------------------------------
#Note that 1000 hr fuels are NOT included in the fuel moisture files

## FUEL MOISTURE - HERBACEOUS FUELS------------------------------------------------------------------------
fmherb.reg <- lm((fmherb~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fmherb.reg)

fmherb.reg.oak <- lm((fmherb.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model
summary(fmherb.reg.oak)

plot(erc.fuel, fmherb, xlim = c(0, max.erc), main="Village Creek Herb. Fuel Moisture")
lines(erc.fuel, predict(fmherb.reg), col="blue")

plot(erc.fuel.oak, fmherb.oak, xlim = c(60, max.erc), main="Oak Knoll Fuel Moisture- High ERC")
lines(erc.fuel.oak, predict(fmherb.reg.oak), col="blue")

pred.fmherb <- predict(fmherb.reg, newdata = newdat) #Predict new herb. fuel moistures with erc.mc1 data
pred.fmherb.oak <- predict(fmherb.reg.oak, newdata=newdat.high)

plot(newdat[,1], pred.fmherb, main="Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted herb. fuel moisture
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fmherb.oak, main="Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted")

#Fuel moistures above 80 are incorrectly predicted!  These data will be replaced.
#Fuel moistures should reach asymptote.
#Replace fuel moisture values where ERC >=80 with 2 based on historical Oak Knoll data
pred.fmherb.oak2 <- ifelse(hightempx>=80, 2, pred.fmherb.oak)
plot(newdat.high[,1], pred.fmherb.oak2) #Plot mc1 ERC by corrected fuel moisture

## FUEL MOISTURE - WOODY FUELS------------------------------------------------------------------------
fmwood.reg <- lm((fmwood~erc.fuel+I(erc.fuel^2))) #Polynomial linear regression model
summary(fmwood.reg)

fmwood.reg.oak <- lm((fmwood.oak~erc.fuel.oak+I(erc.fuel.oak^2))) #Polynomial linear regression model
summary(fmwood.reg.oak)

plot(erc.fuel, fmwood, xlim = c(0, max.erc), main="Village Creek Woody Fuel Moisture")
lines(erc.fuel, predict(fmwood.reg), col="blue")

plot(erc.fuel.oak, fmwood.oak, xlim = c(60, max.erc), main="Oak Knoll Fuel Moisture- High ERC")
lines(erc.fuel.oak, predict(fmwood.reg.oak), col="blue")

pred.fmwood <- predict(fmwood.reg, newdata = newdat) #Predict new woody fuel moistures with erc.mc1 data
plot(newdat[,1], pred.fmwood, xlim = c(60, max.erc), main="Predicted Fuel Moisture - Village Crk",
     xlab="ERC predicted") #Plot mc1 ERC by predicted woody fuel moisture

pred.fmwood.oak <- predict(fmwood.reg.oak, newdata = newdat.high) #Predict new woody fuel moistures with erc.mc1 data
plot(newdat.high[,1], xlim = c(60, max.erc), pred.fmwood.oak, main="Predicted Fuel Moisture for High ERC - Oak Knoll",
     xlab="ERC predicted") #Plot mc1 ERC by predicted woody fuel moisture

#Fuel moistures above 87 are incorrectly predicted!  These data will be replaced.
#Fuel moistures should reach asymptote.
#Replace fuel moisture values where ERC >=87 with 60 based on historical Oak Knoll data
pred.fmwood.oak2 <- ifelse(hightempx>=87, 60, pred.fmwood.oak)
plot(newdat.high[,1], pred.fmwood.oak2) #Plot mc1 ERC by corrected fuel moisture
                        
#------------------------------------------------------------------

#Create dataframe with erc.mc1 and all fuel categories
newfuel <- cbind(erc.mc1, pred.fm1, pred.fm10, pred.fm100, pred.fmherb, pred.fmwood)
newfuel<- as.data.frame(newfuel)
lowfuel <- subset(newfuel, newfuel$erc.mc1 < 60)

highfuel <- as.data.frame(cbind(erc.mc1=hightempx, pred.fm1=pred.fm1.oak, pred.fm10=pred.fm10.oak2,
                                pred.fm100=pred.fm100.oak, pred.fmherb=pred.fmherb.oak2, 
                                pred.fmwood = pred.fmwood.oak2))

##Calculate mean value for each fuel moisture category by ERC value, and then round
fuel.mean1 <- round(aggregate(lowfuel, list(round(lowfuel$erc.mc1, digits=0)), mean), digits=0)
fuel.mean2 <- round(aggregate(highfuel, list(round(highfuel$erc.mc1, digits=0)), mean), digits=0)

##ACTIVATE THE FOLLOWING TWO LINES OF CODE TO VIEW FUEL MOISTURE VALUES BY ERC!!!
fuel.mean <- rbind(fuel.mean1, fuel.mean2) #Create table of fuel moisture values
#write.table (fuel.mean, file="fuelmoisture_hadley.txt", sep=",", quote= F, row.names=FALSE) 

#Set the location to output the files:
## ---> Specify climate model!

path <- "J:\\runfiregen\\fuelmoisture\\HAD\\"

filect <- nrow(fuel.mean)

#Modification Created 11.09.2012
#-----------------------------------------------------
#Create modified fuel moisture files!!!!
#Fuel moisture values for AG-Grass (101, 102, 103, 104 & 107) are HARD-CODED in some cases
#     Added values for Fuel Types 40-69 for non-July files and are hard coded in some cases
#---------------------------------------------------
#Create fuel moisture files, one file per ERC value-- These will apply to all days EXCEPT for exceptions below
#101 - Managed pasture and hay
#102 - High quality savannah
#103 - Structural savannah
#104 - Unmanaged pasture and unmanaged savannah
#107 - Grass seed and grains
#40-69
##7.18.12 Fuel models 101 and 107: 1 hr and 10 hrs fuel moistures increased per request of B. Johnson.
##7.23.12 Fuel models 101 and 107: herbaceous and woody fuels moistures increased per request of A. Ager.
##11.09.12 Fuel models 40-69 added per request of B. Johnson
##11.09.12 Fuel moisture values corrected so base values are always used if higher than modified values
##Note: This coding is not the most efficient, but quickly modified original code

#Insert code to force FM values to base values if preset values are HIGHER Than base values

for (i in 1:filect){
  
  ID <- fuel.mean$erc.mc1[i]
  
  data.now <- paste(paste(("0"), " ",
                          paste(fuel.mean$pred.fm1[i]), " ",
                          paste(fuel.mean$pred.fm10[i]), " ",
                          paste(fuel.mean$pred.fm100[i]), " ",
                          paste(fuel.mean$pred.fmherb[i]), " ",
                          paste(fuel.mean$pred.fmwood[i]), "\n",
                          paste("101"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 12) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >15) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >20) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >150) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("102"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 12) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >15) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >20) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >150) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("103"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 12) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >15) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >20) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >150) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("104"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 12) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >15) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >20) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >150) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("107"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 12) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >15) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >20) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >150) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("40"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("41"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("42"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("43"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("44"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("45"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("46"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("47"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("48"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("49"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("50"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("51"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("52"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("53"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("54"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("55"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("56"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("57"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("58"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("59"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("60"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("61"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("62"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("63"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("64"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("65"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("66"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("67"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("68"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },
                    
                    paste(("69"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 9) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >10) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >15) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >136) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    })
  
  write(data.now, paste(path, paste("FlamMapFuelMoistures", ID,".fms", sep=""), sep=""))
  
} # ends loop


#Assign fuel moisture file to each daily ERC value (data row)
fuel.file <- rep(NA, mc1.obs)

erc.mc1.int <- round(erc.mc1, digits = 0)#Create ERC as integer
fmlabel <- as.data.frame(erc.mc1.int)#Create label for fuel moisture file

for (i in 1: mc1.obs){
  
  fuel.file[i] <- paste("FlamMapFuelMoistures", fmlabel$erc.mc1.int[i], ".fms", sep="")  
} #end loop

#Create fuel moisture files, one file per ERC value-- These will apply to ONLY days July15-Oct1

for (i in 1:filect){
  
  ID <- fuel.mean$erc.mc1[i]
  
  data.now <- paste(paste(("0"), " ",
                          paste(fuel.mean$pred.fm1[i]), " ",
                          paste(fuel.mean$pred.fm10[i]), " ",
                          paste(fuel.mean$pred.fm100[i]), " ",
                          paste(fuel.mean$pred.fmherb[i]), " ",
                          paste(fuel.mean$pred.fmwood[i]), "\n",
                          paste("101"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 12) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >15) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >20) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >150) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n");
                    },          
                    
                    paste(("102"), " ", paste(fuel.mean$pred.fm1[i]), " ",
                          paste(fuel.mean$pred.fm10[i]), " ",
                          paste(fuel.mean$pred.fm100[i]), " ",
                          paste(fuel.mean$pred.fmherb[i]), " ",
                          paste(fuel.mean$pred.fmwood[i]), "\n",  #Predicted values
                          paste("103"), " ", paste(fuel.mean$pred.fm1[i]), " ",
                          paste(fuel.mean$pred.fm10[i]), " ",
                          paste(fuel.mean$pred.fm100[i]), " ",
                          paste(fuel.mean$pred.fmherb[i]), " ",
                          paste(fuel.mean$pred.fmwood[i]), "\n",  #Predicted values
                          paste("104"), " ", paste(fuel.mean$pred.fm1[i]), " ",
                          paste(fuel.mean$pred.fm10[i]), " ",
                          paste(fuel.mean$pred.fm100[i]), " ",
                          paste(fuel.mean$pred.fmherb[i]), " ",
                          paste(fuel.mean$pred.fmwood[i]), "\n",  #Predicted values
                          paste("107"), " "),
                    
                    if (fuel.mean$pred.fm1[i] > 12) {
                      paste((fuel.mean$pred.fm1[i]), " ");
                    } else {paste(("12"), " ");
                    },
                    if (fuel.mean$pred.fm10[i] >15) {
                      paste((fuel.mean$pred.fm10[i]), " ");
                    } else {paste(("15"), " " );
                    },
                    if (fuel.mean$pred.fm100[i] >20) {
                      paste((fuel.mean$pred.fm100[i]), " ");
                    } else {paste(("20"), " " );
                    },
                    if (fuel.mean$pred.fmherb[i] >120) {
                      paste((fuel.mean$pred.fmherb[i]), " ");
                    } else {paste(("120"), " " );
                    },
                    if (fuel.mean$pred.fmwood[i] >150) {
                      paste((fuel.mean$pred.fmwood[i]), "\n");
                    } else {paste(("150"), "\n" );
                    })
  
  write(data.now, paste(path, paste("FlamMapFuelMoistures", ID,"JUL-OCT.fms", sep=""), sep=""))
  
} # ends loop

jdaynum <- as.numeric(jday) #Convert julian day to numeric

#Replace the fuel moisture file name with the JULY15-OCT1 modified files only for Julian days 196-274
for (i in 1: mc1.obs){
  
  if ((jdaynum[i] > 195) & (jdaynum[i] < 275)){
    fuel.file[i] <- paste("FlamMapFuelMoistures", fmlabel$erc.mc1.int[i], "JUL-OCT.fms", sep="")  
  }
} #end loop


###########################
# GENERATE FIRE LIST FILE ##########################
###########################
#Create table that outputs year, julian day of year,
#  fire probability, burn period, wind speed, azimuth and fuel moisture file 
#  for each day in MC1 input file

## year, prob, julianday ,burn period, wind (mph), azimuth, FMS File, daily ERC
fire.list = data.table (yr=year-2007, prob = pred.fire.mc1, julian=jday, burnper = burnp, windspeed = w.spd,
                         azimuth = azimuth, FMfile = fuel.file, ERC = erc.mc1, size = fire.size)

## ---> Specify climate model!

write.table (fire.list, file="firelist_fmmod_HAD20full-list.txt", sep=",", quote= F, row.names=FALSE) #For Record Keeping

#####################################################
# GENERATE FIRE LIST FILE WITH NO FIRES UNDER 1 HA ##########################
#####################################################
#Create table that outputs year, julian day of year,
#  fire probability, burn period, wind speed, azimuth and fuel moisture file 
#  for each day in MC1 input file after implementing ignition switch

#Create data frame of fileid (chardate), erc, ignition, windspeed, azimuth, burn period and fuel file
fire.list.subset <- as.data.frame(cbind(yr=year-2007, prob=pred.fire.mc1, julian=jday, burnper=burnp, windspeed=w.spd, 
                                        azimuth, FMfile=fuel.file, ERC=erc.mc1, size = fire.size), stringsAsFactors=FALSE)
#Grab only ignitions

lrg.fires.only <- as.data.frame(fire.list.subset[fire.list.subset$size >= 1, ])

## ---> Specify climate model!

write.table (lrg.fires.only, file="firelist_fmmod_HAD20.txt", sep=",", quote= F, row.names=FALSE) #For Envision

################################################
# GENERATE FIRE LIST FILE WITH IGNITION SWITCH ##########################
################################################
#Create table that outputs year, julian day of year,
#  fire probability, burn period, wind speed, azimuth and fuel moisture file 
#  for each day in MC1 input file after implementing ignition swtich

#Determine ignition! ##Ignition is determined in Envision- 
#RUN ONLY WHEN YOU NEED PROBABILITY

out.obs <- nrow(lrg.fires.only)

rand.prob <- (runif(out.obs, 0, 1)) #Generate random number
lrg.fires.only$ignition <- ifelse ((lrg.fires.only$prob >= rand.prob), 1, 0)

#Grab only ignitions
fires.only <- as.data.frame(lrg.fires.only[lrg.fires.only$ignition==1, ])
fires.only2 <- fires.only[,-10] 

## ---> Specify climate model!

write.table (fires.only2, file="firelist_fmmod_HAD20_ig.txt", sep=",", quote= F, row.names=FALSE) #For Comparison with Envision

max(erc.mc1) #check max erc
(nrow(fires.only))/98 #check rough estimate of fires/yr
max(pred.fire.mc1) #check ignition probability
max(fire.size) #check max fire size

